Dynamics of electronic transport in a semiconductor superlattice with a shunting side 

layer 
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We study a model describing electronic transport in a weakly-coupled semiconductor superlattice 
with a shunting side layer. Key parameters include the lateral size of the superlattice, the con- 
nectivity between the quantum wells of the superlattice and the shunt layer, and the conduction 
properties of the shunt layer. For a superlattice with small lateral extent and high quality shunt, 
static electric field domains are suppressed and a spatiaUy-uniform field configuration is predicted 
to be stable, a result that may be useful for proposed devices such as a superlattice-based TeraHertz 
(THz) oscillators. As the lateral size of the superlattice increases, the uniform field configuration 
loses its stability to either static or dynamic field domains, regardless of shunt properties. A lower 
quality shunt generally leads to regular and chaotic current oscillations and complex spatio-temporal 
dynamics in the field profile. Bifurcations separating static and dynamic behaviors are characterized 
and found to be dependent on the shunt properties. 

PACS numbers; 



I. INTRODUCTION 

Theoretical work by Esaki and Tsu^ in 1970 was the 
first to propose a Bloch oscillator based on a superlat- 
tice (SL) structure. In that paper, they derived current- 
voltage (I-V) characteristics of a SL which showed neg- 
ative difi^erential conductivity (NDC) associated with 
Bloch oscillations^i^ of the miniband electrons under a 
DC bias. However, direct observation of Bloch oscil- 
lations is difficult due to decoherence caused by elec- 
tron scattering. In other important early work, Kti- 
torov, Simin and Sindalovskii^ predicted a negative high- 
frequency differential conductivity and associated ampli- 
fication of high frequency signals thereby suggesting an 
alternative means of THz oscillation. This dynamic con- 
ductivity remains negative up to the Bloch frequency lob 
and reaches a resonance minimum at a frequency closely 
below ub, suggesting that the SL may serve as an active 
medium for THz radiation. 

However, no such devices have been realized to date 
more than three decades after their proposal because the 
NDC causes space-charge instability. Although Bloch os- 
cillations have been observed experimentally in undoped 
SLs^ by studying optical dephasing of Wannier-Stark 
ladder— excitations using degenerate four-wave mixing^i^ 
the gain achievable is too small to build electrically ac- 
tive Bloch oscillators. For high current densities, the 
space-charge instability causes moving charge accumula- 
tion layers (CALs) and charge depletion layers (CDLs) 
and thus the SL exhibits oscillations similar to the Gunn 
effect.^ While devices based on these oscillations may op- 
erate in the microwave range, they do not extend to the 
THz regionji^ 

The lack of suitable THz radiation sources and detec- 
tors hampers the technological exploitation of the fre- 
quency regime spanning from 300 GHz to 10 THz. Quan- 




FIG. 1: Schematic of the shunted SL. The growth direction 
is along the z direction and the quantum wells are parallel to 
the X direction. The SL is located at a; > and the shunt is 
at a; < 0. The thick line on the right is the potential energy 
of an electron in the conduction band of the SL. 



tum cascade laser devices have been shown to oper- 
ate in the THz range for temperatures up to 164 Kjii 
On the other hand, if superlattice-based Bloch oscilla- 
tors could be successfully realized they might be ex- 
pected to have certain advantages relative to the quan- 
tum cascade structuresji^ Recently, rapid progress in 
THz technologyii including biomedical sensing, three- 
dimensional imaging and chemical agent detection has 
attracted renewed attention to Bloch oscillators. Some 
structures have been proposed to stabilize the field in the 
SL against NDC-related instabilities. One scheme theo- 
retically proposed by Hyart et al}^ is the dc-ac-driven 
SL which requires the presence of an initial THz pump. 
The SL is biased in the NDC region under a DC electric 
field, initially superposed with an AC pump electric field 
which stabilizes the field distributioUfi^ Then the initial 
pump field can be gradually turned off when THz oscil- 
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FIG. 2: Charge density plots (left column) and current vector 
plots of j±m{x) (right column) for a SL with ~ 20 ^im, 
U = 2.1V and a = 0.04 (^m)"^ at 0.02, 0.2 and 6 ns. Initial 
condition is a CAL at the center of the SL. The shunt is at the 
bottom. The color bar on the left of the first contour plot is 
the scale encoding in units of No used throughout the paper. 



lation has been already established in the SL. Another 
suggestion is to stack a few short SLs, where domains 
are not able to forniji^ These short SLs are separated 
by heavily doped material, and an increase in terahertz 
transmission at dc bias has been observed. 

Yet another scheme is to open a shunting channel par- 
allel to the SL, similar to a method that has been used to 
stabilize tunnel diode circuitst^ii^ Daniel et ali^ used a 
distributed nonlinear circuit model to simulate the elec- 
tric field domain suppression in a SL. They have shown 
that the shunt is able to suppress the voltage inhomo- 
geneity above a critical bias voltage which depends on 
the shunt width, the SL width, and the shunt resistiv- 
ity. However, the circuit model does not include aspects 
of the electronic tunneling transport that appear to play 
an important role in SL behavior. The model possesses 
only a global coupling since the elements are connected 
in series and the I ~ V characteristic of each element is 
fixed. On the other hand, the SL model has a more com- 
plex structure that has both a global coupling due to the 
applied voltage constraint as well as a nearest neighbor 
coupling arising from the varying charge densities that 
dynamically change the local current density vs. field 
{J — F) characteristics. As a result, the nonlinear circuit 
model of Daniel et al. is not able to exhibit connected 
field domains or current self-oscillations that are observed 
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FIG. 3: The steady state for a SL with = 20 ^im, U = 
2.1 V and cr = 0.04 (fim)"^ at a; = L^: (a) field profile 
(solid line) and charge density (dashed line), (b) The solid 
dots indicate the actual current operation points on the local 
vertical current field characteristics j||m^m+i(F, n™, Wm+i). 



in SL structures both theoretically and experimentally.— 
In similar work by Fell et al.^^ a side layer is grown 
on the cleaved edge of a lightly doped GaAs/AlGaAs SL, 
such that a 2D electron gas is formed at the interface 
between the SL and the side layer. The lightly doped SL 
serves two purposes: (i) to provide a modulated potential 
for the 2D electron gas at the interface so that under this 
periodic potential, the electron gas becomes a surface SL 
with one lateral dimension; (ii) to provide a uniform field 
to this surface SL since a lightly doped SL can maintain a 
uniform field under external bias. While the suppression 
of field instabilities has been reported in this type of SL, 
it is still not clear whether this lateral structure will be 
useful as a THz oscillator. 

In this paper, we study an extension of a well- 
established model of electronic transport in weakly- 
coupled superlattices by adding a shunting side layer. 
Our treatment includes the effect of lateral electronic 
(i.e., horizontal) transport within each of the quantum 
well layers. Here, the vertical electron dynamics is 
associated with sequential resonant tunneling between 
weakly-coupled quantum well layers, rather than mini- 
band transport or Wannier-Stark hopping as occurs for 
strongly-coupled SLs. Although the Bloch oscillator gen- 
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FIG. 4: Steady states: (a), (d) Charge density plots, (b), (e) 
current vector plots and (c), (f) field profile (solid line) and 
charge density (dashed line) tX x = for ~ 160 (left 
column) and Lx ~ 640 /im (right column), respectively, with 
U ^2.1V and a = 0.04 {flm)-K 
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FIG. 5: SL current density Jsl/L^ and snapshots of charge 
density distribution for — 0.8 mm, U — 2.1 Y and a = 
0.04 {Q,m)~^ . The times of the snapshots are marked as solid 
circles in the upper panel. 



erally requires a strongly-coupled SL, the weakly-coupled 
SL has similar NDC features in I-V characteristics and 
similar current self-oscillations occur due to recycling of 
electronic fronts J^i^^ 

In the next section, we establish a two-dimensional 
model for describing the current flow and dynamical elec- 
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FIG. 6: Same as Fig.[5l but with = 1.28 mm. 



trie field profile in a shunted SL. In section III, we discuss 
the extremely different time scales involved in this model, 
which are challenges to numerically solving it. In section 
IV, we numerically explore the effect of a high quality 
shunt on the dynamics of SLs as the lateral size of the SL 
is varied, and show that the uniform field configuration 
is stable, provided that the shunt and shunt connection 
have high enough quality and the SL lateral extent is 
not too great. In section V, we choose a laterally narrow 
SL and study the dependence of the SL dynamics on the 
shunt properties. The transition from a stable uniform 
field configuration to static field domains is found to be 
complex and the bifurcations involved in this transition 
are discussed. The Appendix presents details of the nu- 
merical methods employed. 



II. LATERALLY EXTENDED MODEL OF THE 
SUPERLATTICE WITH SHUNT LAYER 

Weakly-coupled semiconductor superlattices have been 
successfully described by the sequential resonant tunnel- 
ing model over the past several years .^^SiSSiiS^iSl However, 
previous works usually consider only the dynamics along 
the growth (vertical) direction of the SL and ignore the 
dynamics in the in-plane (lateral) direction, i.e., treat 
each period as an infinitely large plane with uniform 
charge density. More recently, Amann et al^ developed 
a theoretical framework which describes both lateral and 
vertical electronic dynamics. Here, we extend this frame- 
work to include the effects of a shunting side layer. 

The structure of the shunted SL is shown in Fig. [1] 
Each quantum well forms a slab that is parallel to the 
X — y plane, with cross sectional dimensions and Ly. 
There are N such quantum wells stacked on top of each 
other in the z direction, sandwiched between an emitter 
layer and a collector layer. The shunt layer is located 
between —5x < x < 0, with thickness 5x- The SL period 
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is I = w + d, where w and d are the width of the quantum 
well and width of the barrier, respectively. The external 
voltage is applied in the z direction, across the emitter 
and the collector. 

Inside the SL, the electrons are localized within one 
quantum well due to the relatively thick quantum bar- 
riers. Furthermore, the electrons are assumed to be at 
local equilibrium and the local two-dimensional charge 
density at time t is denoted by nm{x,y,t), where m is 
the well index, x, y are the in-plane coordinates. The 
charge continuity equation in the SL can be written as: 

e hrnix,y,t) = j||,„_i^„ - jl\m^rn+l " V± ' j_Lm, (1) 



where 



^ dx 



dy' 



(2) 



and j\\m-i~*m denotes the three dimensional vertical cur- 
rent in z direction tunneling through each barrier (units: 
[A/m^]) and jx,„ is the lateral two-dimensional current 
density (units: [A/m]). The electron charge is e < 0. The 
y-depcndcncc is ignored and Eq. ([T]) can be rewritten as: 



J||rn^m+1 



dj±m{x) 

dx 



(3) 



The local vertical tunneling current j\\m^rn+i through 
each barrier is described by the sequential resonant tun- 
neling model which has been derived using different 
methods j^°'^^i^^ in this paper, we have used the same 
form as in Refs. [20||2^ . This tunneling current depends 
on the electric field F||„(a:) across the barrier through 
which the tunneling occurs and the electron charge den- 
sities nm-i{x) and nm{x) in the neighboring quantum 
wells of this barrier. Thus, the tunneling current has the 
functional form: 



3\\r 



(x) 

— Jllm— l^m 

[F\\jn{x),n„i-i{x),n„-,{x)]. 



. (4) 

The tunneling current densities through the emitter 
and collector layers are modeled by Ohmic bound- 
ary conditions, that is, j||o^i(a;) — cr-F||o(a;), and 
j\\N-tN+i{x) — <jF\\i^{x)niq /Nb, with contact conduc- 
tivity cr and two-dimensional doping density in each 
well. 

The lateral dynamics is caused by the in-plane current 
]jLm which consists of a drift part and a diffusion part. 
When the y-dependence is ignored, this becomes 



dx 



(5) 



where F±jn{x) is the in-plane component of the electric 
field at x in well m, ^ is the mobility and Dq is the dif- 
fusion coefficient. The generalized Einstein relation^ es- 
tablishes the connection between /i and Dq for arbitrary 
two-dimensional electron densities including the degen- 
erate regime: 



Do{nm) 



-epo(l - exp[-n„i/(pofcsr)]) 



(6) 



with the two-dimensional density of states po = 
m* /{nh^), where m* is the electron effective mass. Here 
we assume that and -Dq are fixed. 

Both the lateral and vertical currents depend on the 
electrical fields which in turn depend on the scalar poten- 
tial (pmix, y). The potential can be solved by the Poisson 
equation 

A0m(a;,y) = (A_l + A||)(/)„(x, y) = (n™ - Nd), 

t€r€o 

(7) 

(8) 



with 



A^(t)mix) = ^0m(a;), 



Aj[0,„(a;) = , (9) 

where and eg are the relative and absolute permittivity, 
respectively. Then the field can be calculated as 



F\\mix,y) 



(l)m+l{x) - <pm{x) 



F±m{x) = - 



I 

dcj)m{x) 

dx 



(10) 



Here we solve the Poisson equation using an approxima- 
tion method assuming that the typical structures in the 
lateral direction vary on a length scale much longer than 
the mean free path of the degenerate electrons.^^ 

The drift-diffusion dynamics of the shunting layer is 
similar to that of the lateral dynamics within each SL 
quantum well. First, we neglect x-dependence in the 
shunt, that is, the shunt is collapsed into a single layer 
along the z-direction. Note also that unlike the SL, which 
possesses an intrinsic discreteness along z direction, the 
shunt is a continuous layer. Therefore, we make a fur- 
ther approximation that the shunt is divided into blocks 
aligned with the periods of the SL and that the charge 
density is locally uniform within each block. This as- 
sumption not only provides the discretization required 
by numerical simulation, but also matches the dynamics 
of the shunt with that of the SL. With these two assump- 
tions, we can write down the continuity equation in the 
TO-th shunt block as follows: 



e niJ (t)-Sx I Ly = Til"'* , 



(11) 



where the superscript (s) denotes the quantities in the 
shunt and the tilde denotes that the quantities are three- 
dimensional, i.e., 



Here, the quantity j^''^ denotes the lateral current that 
flows between the shunt and the SL through their inter- 
face. Then we can write Eq. pT|) in the form: 
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Note that the vertical current in the shunt has a very 
different form than the tunnehng current in the SL. It 
follows a similar dynamics as the in-plane current in the 
SL quantum wells and is related to the three-dimensional 
charge density in the shunt: 



J\\m-1- 



dfi 



dz 



(14) 



Here we assume the mobility /i and the diffusion coeffi- 
cient -Do have the same values as in the SL. 

Next, we examine the lateral current that connects the 
shunt and the quantum well layer within the SL: 



_Lrn 



(15) 



a;=0+ 



In this equation, the boundary should be defined at 
X = 0+ for calculation of both the current and the po- 
tential in the shunt. Since the shunt is assumed to be 
uniform in x direction, defining the above equation at 

— . . (s) 

X = 0^ implies that F±m and Vj_rim are zero which 
would lead to zero boundary current. Another advantage 
of choosing the boundary at a; = 0+ is that the potential 
in the shunt should be equal to the potential in the SL 
close to its boundary, i.e., ipm {x < 0) = 0m (a: = 0+), 
since the potential is continuous everywhere. This rela- 
tion allows us to equate the potential in the shunt with 
that at the inner boundary of the SL. So the potential 
at the boundary of the solution of Eq. ([7]) is just the po- 
tential in the shunt. The fields required to calculate the 
current in Eq. psp can be obtained by 



i^im(0+) = -Vi0™(x) 



c^'^:i+A^)-c^t\x) 



(16) 



The charge density and its normal gradient at the bound- 
ary are 



njn{x = 0) 
V I n,, 



(17) 



x=0+ 



= hm ^^n{Ax)-ni^\ 

Ax^O+ Ax 



SL can cause band bending effects at the interface. Even 
if the shunt is doped to have the same Fermi level as that 
in the SL so that little band bending might be expected, 
there are other issues that impact the connection quality 
between the shunt and the SL, for example, the presence 
of trap states or a thin oxide layer. To quantify the qual- 
ity of the connection between the SL and the shunt, we 
introduce a parameter < a < 1 such that a = 1 corre- 
sponds to a perfect connection and a — corresponds to 
no connection. We modify Eq. (fT5| to be 



a ■ -enriynix = 0)F±rn - DqV j_nr, 



=0+ 



(19) 

The relationship between specific values of parameter a 
and microscopic models of conduction across the shunt- 
SL interface are discussed elsewhere.— 

Similarly, we introduce a separate parameter 6 > 
that allows us to model the effect of having different dop- 
ing density and/or mobility in the shunt vs. SL quantum 
wells. Also recognize that the field in the shunt is al- 



most uniform and n 
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when the conductance in 
is the doping density in 
This leads to the following modification of 



(s) 

the shunt is high, where N^' 
the shunt. 
Eq. (HI, 



■^llm— 1 — >m 



eDn 



dUm 



As) 



(20) 

where bfiNo ^i'-'^N'i^\ Note that b > 1 when the 
doping density in the shunt is greater than that in the 
quantum wells and b is much less than one when the 
shunt is weakly conducting so that only a small fraction 
of the total vertical current flows through it. 

It is also useful to point out that the total current. 



J 



is) 



As) 



m — *m4-l 



y^xA J (cr-eo^'llm + j|lm^m+l) dx, 

(21) 



is the same for each period. To show this, note that the 
Poisson equation can be written as 



V- (Fj 



lerCo 



[Urn - Nd), 



(22) 



Here we also note the possible effects of energy band 
structure of the shunted SL and the doping density in the 
shunt. In the above discussion, the situation has been 
simplified because no band bending is included. How- 
ever, variations in doping densities in the shunt and the 



|m-l 



dx hreo 



(n„ 



N 



D) 



(23) 



Substituting the above equation into Eq. ([3]) yields 
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Then, one integrates both sides of the preceeding equation with respect to x from —Sx to L^- Due to the vanishing 
boundary conditions F±(—Sx) — F±{Lx) — and j±m{~^x) = j±miLx) — 0, the lateral terms in the above equation 
integrate to zero. This yields 

erCo^y F\\mdx + J j\\m^m+idx = ereo-^ J F||,„_idx + y j||„_i^„da;. (25) 



which shows that the total current is independent of the 
well index m. Note that the current through the shunt 
will be the dominating contribution to the total current 
of a SL if the shunt is thick and well-conducting. Even 
a completely disconnected shunt (i.e. a — 0) contributes 
a constant current of Jq*^ = SxCplNdU / {Nl + d) to the 
total current J of a homogeneous SL. Since we are inter- 
ested in effects arising from the interaction between the 
SL and the shunt, we will in the following discuss the 
current dynamics on the basis of the SL current defined 




III. PARAMETERS AND TIME SCALES 

The parameters that we use in the simulation are listed 
in Table L We found that there are very different time 

TABLE L Parameters used for the shunted SL. 
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FIG. 7: Same as Fig. [3 but with = 2.56 mm. 



scales in this complex structure which requires an implicit 
method of numerical iteration. The first time scale Tf, is 
the dielectric relaxation time in the bulk material both 
in the shunt and in each quantum well in the SL. It is 
determined by the doping density. We know that the 
conductivity g is proportional to the charge density 

g « enNo/l ^ 1.6xl0"^^xl0xl02^(f7m)"^ - 10^(f7m)"^ 

(26) 

So the dielectric relaxation time in the shunt layer and 
within each quantum well is approximated as 



n 



9 



0.1 X 10"^ 
10^ 



(s) ^ lO-i^(s) 



(27) 



which is relatively fast due to the high conductivity. This 
is the time it takes for a fluctuation in the charge den- 
sity to be neutralized within either the shunt layer or 
quantum wells. 

The second time scale tj is the one in the vertical dy- 
namics. According to the sequential resonant tunneling 
model, the vertical current is to the order of 10~* (A/m^) 
and the positive differential conductivity gt is of order 
0.1 (17 m)~^. Thus, Tt = erCo/gt ^ 10^^ s, a much larger 
time scale than Tf,. Moreover, from numerous previous 



works, we also know that the behavior of the electrons 
in the vertical direction is not simply dielectric relax- 
ation. More complex phenomena, such as current self- 
oscillation, or injected dipole relocation due to switching, 
have much longer time scales ranging up to microseconds. 
The time scale Tt sets a lower limit of the time scales for 
these nonlinear processes. 

Another important time scale is the time that it 
takes to carry away or supply the electrons in the SL 
through the shunt. Because the vertical processes are 
relatively slow, if the shunt has good connection and 
high conductance, the electrons will move laterally, pass 
through the intersection between the quantum well and 
the shunt, and drift away through the shunt. This time 
scale Ti is considerably larger than Tf, since the electrons 
have to move into the shunt first. Later we will see that it 
takes 1 ns to deplete a full CAL in a small SL. The pres- 
ence of extremely different time scales means that the 
numerical integration is a stiff problem and this suggests 
the use of an implicit method. The numerical procedure 
is described in the Appendix. 
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FIG. 8: Same as Fig. [H but with = 5.12 mm. 



DEPENDENCE OF SHUNTING DYNAMICS 
ON THE LATERAL SIZE OF THE 
SUPERLATTICE 



In this section, we discuss the effects of the lateral size 
Lx of the SL with a high quality shunting layer, i.e., 
a — b = 1. The shunting layer has a width 6x such that 
varying 6x does not affect the dynamics in the shunt. 
This is numerically confirmed even for the chaotic case 
that we will discuss below, where a 80 nm shunting layer 
has the same effect as a 8 mm one. This is because Tb is 
much smaller than and the electrons entering the shunt 
are carried away so fast that a change in the shunt con- 
ductance does not change . We will study the SLs with 
a relatively high contact conductivity cr — 0.04 (fim)^^. 
At this value of cr, without a shunt, the SL has a static 
high field domain near the emitter and a static low field 
domain near the collector separated by a static charge ac- 
cumulation layer (CAL). Due to the high quality shunt 
the total current is dominated by the contribution of the 
current through the shunt. As discussed at the end of 
Section HH we will therefore consider the SL current Jsl- 
Also, since we are varying L^, we scale current to current 
density. 
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FIG. 9: (a) Bifurcation diagram for a — 0.04 (flm)^^, 
Lx = 20 fim, b = 1.00. Dashed curve shows the approxi- 
mate boundary of the oscillatory region and location of stud- 
ied bifurcation points A and B. (b) Bifurcation scenario at 
A for a = 1.00 x 10~^: amplitude vs. voltage (main figure) 
and frequency vs. voltage (inset). Points Au and Ai denote 
the endpoints of the upper and the lower branches, respec- 
tively, (c) Bifurcation scenario at B for a = 1.00 x 10~^: 
scaling of frequency vs. voltage (double logarithmic plot); 
UcritB = 2.30441 V. 



A. High quality shunting layer with small 

Figure [2] shows charge and current density plots for a 
relatively narrow SL with lateral extent — 20 fiia. 
The initial state is prepared as a charge configuration for 
the SL without shunt at total applied voltage U — 2.1 V 
and shows a static charge accumulation layer at the 20th 
period. After an interval of about 1 ns, the space charge 
configuration is almost uniform. The in-plane current is 
plotted as a vector field and shows the electrons in the 
CAL move in the lateral direction (the opposite direction 



of the current) into the shunt. We can see that when 
the system reaches steady state, the net charge is almost 
neutral, i.e., n = Njj, everywhere in the SL and the 
shunt. There are still some small lateral current flows at 
the first and the last period. 

If we take a close look at the steady state, we find 
that there is a small CAL at the first period and a CDL 
nearby (Fig. [3ja)). The situation is almost inverted near 
the collector. To better understand this, we focus on the 
operation points near the emitter shown in Fig. ^h) at 
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FIG. 10: (a) Charge density distribution evolving in time, 
(b) a snapshot of field profile (solid line) and charge density 
profile (dashed line) at t = 8 ns and (c) SL current Jsl on 
the upper branch of Fig. [9l3. Parameters: a — 1.00 x 10~^, 
U = 1.46 V, a = 0.04 (J7m)-\ 



X = 20 fim. In this case, the field is almost uniform in 
the SL and each period is biased in the NDC region. The 
field across the first barrier between the emitter and the 
first well will also have this same value in the absence 
of charge accumulation in the first well. This causes a 
vertical current from the emitter to the first period (thin 
solid line in Fig. [3fb)) which is much larger than the 
vertical current in the corresponding NDC region of the 
SL. Close to the shunt this extra current will give rise to 
a lateral current which will quickly reach the shunt and 
is carried away by the shunt. A little further away from 
the shunt where the lateral current is not sufficient to 
completely neutralize this extra current, a small CAL is 
formed in the first well which lowers the electric field and 
therefore the current across the first barrier. At the same 
time, the electric field in the second barrier is pushed 
above the uniform field, causing a very small CDL next 
to the CAL. Similar arguments can be applied to the 
collector to explain the appearance of a small CDL in 
the last quantum well. The overall effect is that a nearly 
uniform vertical electric field configuration is stabilized 
for these conditions. 



B. High quality shunting layer with large Lx 
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FIG. 11: (a) SL current Jsl vs. time for U = 1.486 V (the 
inset shows an enlargement), (b), (c) charge density distribu- 
tions evolving for two different time intervals for voltage near 
bifurcation point Ai. (d) SL current for U — 1.481 V. (e) 
The rate A of exponential decay A < (or increase A > 0) of 
oscillation amplitude versus applied voltage U. Parameters: 
a = 1.00 X 10"^ a = 0.04 (nm)-^ 



As the lateral size of the SL becomes larger, the 
CAL and CDL near the emitter become more prominent 
(cf. Fig. [Ha)-(c), Lx = 160 /xm) since with increasing 
distance to the shunt the lateral current becomes less 
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FIG. 12: (a) SL current Jsl (the inset shows an enlarge- 
ment), and (b) charge density distribution evolving in time 
for U = 1.49985 V with a = 1.00 x 10-^ a = 0.04 {flmyK 
(c) Power spectrum data for oscillations ai U — 1.4997 V and 
U = 1.49985 V. (d) The time T for which the system exhibits 
transient oscillations versus the applied voltage U — Ua^ ■ 
Ua,, = 1.499791 V. 



efficient at carrying away the excess current from the 
emitter to the shunt. 

For wider SL (cf. Fig. |4i;d)-(f), = 640 /.tm), the 
field closer to the shunt is more uniform and the CAL 
is still attached to the emitter. However, away from the 
shunt, the CAL detaches from the emitter and locates 
itself in the first few periods and the nonuniform field 
region becomes larger. This behavior is due to the lateral 
current being insufficient to carry away the extra current 
from the emitter. Thus, the CAL grows bigger and tends 
to move toward the collector. With the center of the 
CAL located in different wells at different x positions, the 
lateral gradients can be increased and a sufficient lateral 
current can be sustained. The field profile at a; = 640 /im 
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FIG. 13: Bifurcation scenario at point B of Fig. 9: (a), (c) 
charge density distributions vs. time, (b), (d) SL current Jsl 
for U — 2.3 V (upper panel) and U = 2.304 V (lower panel), 
respectively, with a = 1.00 x 10"^ a = 0.04 (!:2m)"\ 
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FIG. 14: Charge density distribution vs. time for a = 1.00 x 
10"^ a = 0.04 (nm)-i at (a) U = 1.2 V, (b) f/ = 1 V, and 
(c) U = 0.5 V. (d) SL current at (/ = 0.5 V. 



is plotted in Fig. Field domains are forming as 

the field is low to the left of the CAL and high to the 
right of the depletion region. In this case, the upstream 
CAL (closer to the emitter, at the left bottom corner 
of Fig. IHd)) and the downstream CAL (closer to the 
collector, the wider one in Fig. [Hd)) are still connected 
and this is a time-independent steady state. 

In the above case, the lateral size of the SL is just 
below a characteristic value for which the steady state 
loses stability to oscillatory behavior. Figure [5] {L^ ~ 
800 /im) shows the simulations of a slightly wider SL than 
considered above. The large downstream CAL still stays 
in that position. However, due to the large size of the 
SL, the lateral current is not able to sustain a connected 
stable CAL. The small upstream CAL touches and breaks 
off from the downstream CAL periodically. There is a 
small amplitude oscillation in the total current which is 
shown in the top panel. 

For an even wider SL (Fig. [6] with — 1.28 mm), 
the upstream and downstream CALs are mostly discon- 
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FIG. 15: (a) Bifurcation diagram for a — 0.016 (f2m) ^. 
Charge density distribution vs. time near bifurcation point 
C, a = 1.00 X lO-"* at (b) U = 1.6 V, (c) U = 2.1 V, (d) 
U = 2.7 V, and (e) U = 3.0 V. 



nected. The upstream CAL extends laterally into the SL 
and moves toward the downstream CAL, (at time 1.969 
ms). For certain times during the dynamical evolution 
(not shown in Fig. (5]), the upstream CAL breaks off from 
the emitter and reaches and merges with the downstream 
CAL. Mostly, there is a depletion region forming between 
the upstream and downstream CALs (2.211 ms). This 
depletion region grow and dies away as it merges with 
the upper CAL (1.969 ms). For certain times, it grows 
into a full CDL extending throughout the structure and, 
in this case, the upstream CAL also grows into a full CAL 
(1.755 ms). Then all three fronts move downstream. The 
old downstream CAL and the CDL quickly disappear and 
the new CAL formed by the upstream CAL replaces the 
old CAL and stays at its position. Although these be- 
haviors are quite complicated, they are still periodic and 
during each period, the upstream and downstream CALs 
merge several times. 

However, for an extremely wide SL (Fig. ^ = 
2.56 mm), the behavior is apparently chaotic. The effect 
of the shunt is to cause a CAL attached to the emitter 
near the shunt. For large values of x, the shunt has less 
effect and this CAL detaches from the emitter, tends to 
move downstream to the collector and thus extends to- 
ward the downstream CAL. Due to the large lateral size 
of the SL, the impact of the shunt layer becomes very 



weak on the opposite side of the SL. Thus, the down- 
stream CAL is located very close to the 20th period where 
it would be in the absence of a shunting layer. The merg- 
ing of the CALs described in last paragraph also appears 
here except that the merging events are now difficult to 
predict and manifestly not periodic. Figure [5] shows the 
behavior of a SL with = 5.12 mm. It should be noted 
that real SL samples rarely have such a large size. In this 
case, the unstable dynamics only occurs in the portion 
of the SL closest to the shunt. In the portion of the SL 
away from the shunt, a CAL is located at the 20th well, 
where the shunt has no apparent influence. Over time, 
the lateral extension of this CAL changes. When a large 
CDL collides with it at 5.696 ms, the static CAL shrinks 
to a small size, causing a large dip in the current trace. 
The presence of such charge tripole configurations'^* of 
one CDL and two CALs has already been shown to be 
associated with chaotic behavior in one-dimensional SL 
models without lateral dynamics. 

To summarize, we are able to identify three character- 
istic length scales in the x direction. The shortest one 
is the decay length (of order 10 /xm) at which the 
charge density in the first quantum well increases from 
No at the SL-shunt interface to its maximum value (cf. 
Fig. [Ha)-(c)). The next length scale (of order 200 ^m) 
is the range above which the vertical field configuration 
loses uniformity and static field domains start to form 
(cf. Fig. |D(d)-(f)). The longest length scale (of order 
700 /im) is the width of the SL above which the steady 
state loses stability to oscillatory behavior. This implies 
that lateral uniformity in the electric field distribution 
can be expected when is smaller than the intermedi- 
ate characteristic length scale. The shortest decay length 
Lx can be estimated by noting that the extra current 
coming from the emitter must be directed to the shunt 
by the negative gradient of the lateral current J±, i.e., 
^'^g^^^ = J||o^i(a;) — J|| 1^^2(2;) < 0. Then there is ap- 
proximately a decay length L^, at which the quantities 
such as J±{x), n(x) and Fx(x) approach asymptotic val- 
ues exponentially. Calculation shows that is of order 
10 /im for the parameters used in Table 1^ in agreement 
with our numerical results. 



V. DEPENDENCE OF DYNAMICAL 
BEHAVIOR ON THE SHUNT PROPERTIES 

In the previous section, we have seen that the width 
of the SL determines the lateral dynamics of electronic 
transport and that the shunt can stabilize a nearly uni- 
form field configuration in sufficiently narrow SLs. Now 
we investigate the effects of the shunt properties on a 
small SL with width of 20 /im where the lateral field 
and electron density profiles are almost uniform. Since 
the charge density is almost uniform laterally, we modify 
the model such that the SL is collapsed to one point in 
X direction. This modification significantly reduces the 
complexity of the simulation. We first study the effects 



11 



of connectivity parameter a on a SL with conductivity 
cr = 0.04(rini)~^ chosen as in the previous section. Then 
we study the effects of a on a SL with lower contact 
conductivity a ~ 0.016(ilm)^^, which corresponds to 
moving fronts and current self-osciUations in unshunted 
SLs,— and briefly discuss the effects of shunt conductiv- 
ity parameter b and width S^. Since Lx is fixed, we plot 
the unsealed SL current Jsl- 



A. Dynamical behavior vs. connectivity parameter 
a for large contact conductivity 

Figure [9ja) shows a bifurcation diagram using as the 
bifurcation parameters the connectivity parameter a and 
the voltage U for a = 0.04(rim)~^. There is a bounded 
region where the system exhibits periodic or chaotic os- 
cillations, shown as the region enclosed by dashed lines 
in Fig. [ni^a). The value of the connectivity parameter a 
of the oscillatory region ranges from about 6 x 10"'^ to 
7 X 10~^. In real samples, such a weak connection be- 
tween the SL and the shunt could be associated with a 
potential barrier formed between the SL and the shunt 
due to band bending or an oxide layer. 

For a > 6 X 10"'^, the charge density in the SL is 
almost uniform except for a small CDL near the emitter, 
the same situation shown in Fig. [2l With the increase 
of voltage, this CDL becomes more prominent and there 
is an CAL in the first period. However, this CAL never 
detaches from the emitter for any value of voltage when 
a ^ 6x10"'^. This is reasonable because fora>6xl0~"^, 
the connection is strong enough that the shunt is able to 
maintain the field in the SL almost uniform. 

Another stable region is a < 7 x 10^^. In this region, 
a static CAL is formed in the SL and located close to 
the position where it is expected when there is no shunt. 
This is also easy to understand because the connection 
is so weak that the shunt has almost no influence on the 
SL. 

Between these two values of a, we have a transition 
region where oscillations occur for certain ranges of ap- 
plied voltage. Here, the bifurcation scenarios by varying 
voltage are investigated for two sets of values {A and B) 
of the control parameters. 

The bifurcation for point A occurs at a = 1.00 x 10^'^ 
and U = 1.485 V (Fig.[Wb).[Tmfn|. Inside the oscillatory 
region (approximately U < 1.485 V), the charge density 
distribution in the SL oscillates and the oscillation only 
involves part of the SL (cf. Fig. [rU]) . There is a static 
CAL near the emitter but this is clearly detached from 
the emitter. The oscillation occurs in the wide region 
to the right of the CAL in the form of moving charge 
dipoles (CALs and CDLs), cf. Fig. [TUf a). However, at 
any given time, there are three to four pairs of dipoles 
present. From Fig. [TOT b). we can see that the charge 
densities have large amplitude fluctuations along the z 
direction and the higher frequency component of the cur- 
rent oscillation is due to the movement of these dipoles 



(Fig. UnUc)). This higher frequency /i is nine times the 
lower one /2 at which the collector receives the moving 
dipoles. Here we observe the coexistence of static CAL 
and steady moving fronts. 

The bifurcation scenario of A is illustrated by Fig.[9l^b), 
where the amplitude of the current oscillation is plotted 
versus the applied voltage. There is a bistability region 
between U « 1.485 V and 1.50 V, where the system either 
oscillates (upper branch) or is in a steady state (lower 
branch) . 

The bifurcation at point Ai at the end of the lower 
branch is studied in Fig. [TT] When the system starts 
from a uniform configuration ai U ~ 1.486 V, shown in 
Fig. fTlTa'). (b), (c), it first oscillates similar to the full 
oscillation in Fig. (TUl except that the CALs and CDLs 
are much smaller in Fig. [TlTb). The oscillation gradu- 
ally decays to a steady state where there is only a single 
stable CAL and no charge fronts to its right, as shown 
in Fig. [TlTc). The amplitude of the current oscillation is 
quite small and decays to zero. The well-to-well hopping 
of the small charge fronts does not have an appreciable 
effect on the current oscillation form as found for the ma- 
ture fronts in Fig. fnTT cV Instead, the shape of the current 
oscillation is smooth and sinusoidal and possesses a well- 
defined frequency. After a transient interval, the ampli- 
tude A{t) of the current oscillation decays exponentially, 
i.e., A{t) = A{to)expXt and the rate A can be deter- 
mined by fitting. It also should be mentioned that the 
initial state corresponding to the uniform field configura- 
tion falls into the basin of attraction of the upper oscilla- 
tory branch for U ^ 1.486 V. Hence, to obtain A for the 
lower branch, we start the system from the steady state 
oiU ~ 1.486 V. This initial state is used for all the points 
of the lower branch. In the case oi U = 1.481 V, shown 
in Fig. Illf d). the amplitude of the current oscillation in- 
creases exponentially at first and after passing a certain 
threshold value, quickly evolves into the large oscillations 
of the upper branch. The inset shows the transition re- 
gion and indicates that the small charge fronts grow into 
mature ones. The rate A can also be fitted and now it is 
positive. The resulting A versus U is plotted in Fig.fTlTe). 
showing a linear scaling. This clearly indicates that the 
bifurcation at Ai is a subcritical Hopf bifurcation. Su- 
percritical Hopf bifurcations in different SL models have 
been found by Patra et a/.^° and by Hizanidis et al^i at 
low contact conductivity with no shunt. Here we can also 
see that the time scales have the following relationship: 

n < i//i,2 < i/A. 

It is likely that the bifurcation scenario at A^ in 
Fig. m^b) is a saddle-node bifurcation which is proba- 
bly caused by the collision of the stable limit cycle and 
the unstable limit cycle that arises from the subcritical 
Hopf bifurcation at Ai. In Fig. [T^c), the power spec- 
trum of the limit cycle {U — 1.4997 V) and the power 
spectrum of the transient oscillation a,t U = 1.49985 V - 
which exceeds the saddle-node bifurcation value Ua^ - are 
almost identical. This rules out a subcritical torus bifur- 
cation. Then we start the system from a configuration 
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corresponding to the steady oscillation at U = 1.46 V, 
but for voltages just above Ua^ where there are no limit 
cycle states, so it eventually reaches the lower branch. 
Figure [Ufa) shows this process at U = 1.49985 V. Af- 
ter a short time interval of about 1 /is, the oscillation 
amplitude A{t) enters a regime of transient oscillations 
and after a relatively long time T, it suddenly exits this 
region and reaches a steady state. This process looks like 
a reverse process of Fig. irTTd) . Figure [T2\ h) shows the 
decay of the CALs and CDLs. If we choose the critical 
value to be 1.499791 V, then the slope in Fig. \W[d) is 
-0.5. This means that T oc , ^ consistent with a 

system that undergoes a saddle-node bifurcation of limit 
cycles. '^^ 

The bifurcation at point S is at a = 1.00 x 10^^ 
(Fig. [13]). For U < 2.305 V, the system oscillates. At 
first, there is a single CAL in the SL and a dipole is in- 
jected from the emitter. The CAL and dipole all move 
into the SL. The leading CDL moves about twice as fast 
as the two CALs^^ and when it catches up with the orig- 
inal CAL, they annihilate. The CAL of the dipole con- 
tinues to move forward until it reaches the position of 
the original CAL and stays there for a certain period of 
time, waiting for another round of dipole injection. Such 
a bifurcation of a stationary domain state has been re- 
ported before by Hizanidis et al.^^ for a one-dimensional 
superlattice model without shunt at higher contact con- 
ductivity. The time needed for a dipole to be injected is 
called the activation time and the time needed to return 
from the excited state to the fixed point is called the 
excursion time.^ As the applied voltage U approaches 
the boundary, the activation time becomes longer and 
longer. Taking the critical value UcritB of voltage to be 
2.30441 V and plotting the frequency of oscillation versus 
U — UcritB 1 we find the frequency obeys the square-root 
law which is the characteristic scaling law for the saddle- 
node infinite period bifurcation or SNIPERf2i which is a 
global bifurcation of a limit cycle. 

Inside the oscillatory region in the a — U parameter 
space (Fig. UKa)), we also find regimes of chaos. We still 
use a = 1.00 x 10~^. As the voltage U decreases inside 
the oscillatory region, the oscillation shown in Fig. [TUl 
involves a larger part of the SL and the CAL near the 
emitter becomes less and less prominent until these mov- 
ing dipoles cover almost all the SL shown Fig. fWal.fb) 
at U = 1.2 V and 1.0 V. Further decrease of the volt- 
age causes the disappearance of the static CAL and the 
dipoles either annihilate inside the SL or reach and disap- 
pear at the collector, shown in Fig. [TW c'l for U = 0.5 V. 
Similar chaotic behavior has also been found in SLs with- 
out a shunt 1^ These complicated and apparently chaotic 
oscillations are found at many points in the oscillatory 
regime of Fig. Ufa). 

In the regime of the stable states between a « 6 x lO^'^ 
and 7 x 10~^ (cf. the right hand region of Fig.[9l^a)), the 
SL usually has a static CAL either inside the SL (for low 
a) or attached to the emitter (for high a) and there is a 
small static CDL to the right of this CAL. This means 



that the overall field profile is nearly uniform for larger a 
(> 6x 10^"^), but static field domains form as a decreases. 



B. Dynamical behavior for small contact 
conductivity 

The bifurcation scenario for lower contact conductivity 
a is simpler than for the high cr case. Figure [15] shows the 
bifurcation diagram for a = 0.016 (fim"^). This value of 
a corresponds to current self-oscillation in the SL when 
there is no shunt.— The parameter space is again divided 
into an oscillatory regime and a stable stationary regime. 
The oscillatory regime starts at about the same value of 
a as the high a case, i.e., a « 6 x lO^'^. However, this 
oscillatory region does not have a lower bound. This is 
because without the shunt the SL still exhibits oscilla- 
tions. 

The different behaviors at a = 1.00 x 10^"* are shown 
in Fig. [15] As the voltage is deep inside the oscillatory 
region, dipoles are periodically injected into the SL and 
travel through the entire SL (Fig. [TsFb)). As voltage 
increases, the distance that the dipoles travel becomes 
shorter and the CAL and CDL annihilate near the emit- 
ter (Fig. [TST c^V Similar behaviors have been found in a 
SL model without shunt.— As the voltage approaches the 
boundary, the CDL becomes less and less prominent and 
the length that the CALs travel becomes even shorter. 
After the voltage crosses the boundary, the CALs be- 
comes static. The bifurcation scenario is similar to point 
A, described in the previous section, where there is bista- 
bility between oscillatory and steady states. The bifurca- 
tion scenarios at other points on the right hand boundary 
of the oscillatory region appear to be similar to that at 
point C. 



C. Dynamical behavior vs. shunt conductivity 
parameter b 

The above discussion focuses only on varying the con- 
nectivity parameter a with a shunt of high conductance. 
It is also possible to change other parameters of the shunt, 
such as the conductivity parameter b in the shunt. A bi- 
furcation diagram can be plotted for b versus U with 
fixed a = 1.00 and Sx = 200 nm, and it is similar to 
that shown in Fig. [9] with an oscillatory regime between 
w 4.5 X 10^* and 4.5 x 10"''. Another possible control 
parameter is the width of the shunt. Simulation shows 
that only when the width of the shunt is narrower than 
about 1 X 10"* nm, which is unrealistically small, the 
SL starts to have oscillation. The oscillatory region for 
a in Fig. [5] and [15] is almost not affected when b and 6x 
are above certain values so that the current between the 
shunt and the SL can always be supported by the shunt. 
In reality, Sx and b should be kept as low as possible to 
reduce the power dissipated in the shunt and minimize 
heat production. 
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VI. CONCLUSIONS 

We have theoretically studied the effect of a shunting 
side layer parallel to a semiconductor superlattice, and 
find that such a structure can have an almost uniform 
electric field over the entire structure even when biased 
in the negative differential conductivity (NDC) region. 
However, even for a shunt with high conductivity and 
strong connection to the SL, the field in the SL can be 
stabilized only for structures with relatively small lateral 
extent. As the lateral size Lx becomes larger, the lateral 
current in the quantum well loses the ability to deplete 
the extra current coming from the emitter and the field 
becomes nonuniform. For a sufficiently thin SL whose 
lateral dynamics is uniform, the connection between the 
shunt and the SL and the conductivity of the shunt deter- 
mines the dynamics in the SL. We have also established 
the bifurcation diagrams for SLs for different values of 
the shunt parameters and identified the presence of both 
local (Hopf) and global (SNIPER) bifurcations. 

Although the microscopic nature of electronic trans- 
port in weakly-coupled SLs is different than for strongly- 
coupled SLs, the NDC property is known to produce sim- 
ilar dynamics in both types of structures when they are 
not shunted. Thus, it seems plausible that for suitable 
shunt connectivity and SL lateral width that a stronly- 
coupled SL might also be stabilized with a shunting side 
layer. This could enable the realization of a SL-based 
THz oscillator. 



lows: after discretization of the space, the quantities of 
potential and charge density are placed on the grid. The 
fields and currents (also the charge density that is needed 
to calculate the currents) are placed on a staggered grid. 
Knowing the charge density distribution, the potential is 
determined by the Poisson equation using a method de- 
scribed in Ref. '25'. After that, the currents to each grid 
point are calculated from the electric fields which are im- 
mediately obtained from the potential (cf. Eqs. (|10p and 
(fT6| ). Then the charge densities are iterated one step 
forward in time as 



en' = en + diJ(n'), 



(A.l) 



where n = (nn, ni2, n2i, '^■22, •■•)"^ the vector whose 
components are the charge densities on each grid point. 
The first subscript denotes the SL period number and 
the second one is the grid point index in the x direc- 
tion. The vector current J is the total current fiow into 
or out of each grid point, n' is the new charge den- 
sity configuration after time step dt. Since we are using 
the implicit method, J must depend on the future charge 
density configuration instead of the old one. We linearize 
the equations: 



n' = n + dt 



J n + — 
on 



(n'-n) 



(A.2) 



where dJ/dn is the Jacobian matrix. Rearranging this 
equation yields: 
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APPENDIX: NUMERICAL METHOD 

In order to implement the implicit method, the dynam- 
ical variables, i.e. the electron densities nm{x), should 
be computed from the system Eqs. ([3]) and (fT3|) . How- 
ever, nm{x) are deeply buried in these equations, where 
the currents depend on the field that relates to nm{x) 
by solving Poisson equation Eq. ([7]). So instead of solv- 
ing for nm{x) directly, we use the semi- implicit Euler 
method and numerically calculate the Jacobian matrix 
that is needed for this method. The procedure is as fol- 



n' = n + dt 



1-dt 



d3_ 
dn 



J(n) (A.3) 



We mentioned that the currents do not depend on the 
charge densities explicitly. So to calculate the Jacobian 
matrix, we first calculate J(n), then slightly change the 
charge density at one grid point to nij+6nij and calculate 
the currents J' based on this charge configuration. Then 
one row of the Jacobian matrix is immediately obtained 
by (J'-J(n))/(5n„. 

To solve Eq. (|A.3p . we do not invert the matrix. In- 
stead, we write it as: 



diJ(n) = 



, d3 
on 



(n'-n). 



(A.4) 



Then we solve this set of linear equations by Gauss elim- 
ination. 
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